arXiv:1507.06426v2 [cond-mat.quant-gas] 20 Dec 2015 


Properties of the one-dimensional Bose-Hubbard model from a high-order 

perturbative expansion 

Bogdan Damski and Jakub Zakrzewski 
Instytut Fizyki imienia Mariana Smoluchowskiego, 

Uniwersytet Jagiellonski, ulica Lojasiewicza 11, 30-348 Krakow, Poland 

We employ a high-order perturbative expansion to characterize the ground state of the Mott 
phase of the one-dimensional Bose-Hubbard model. We compute for different integer filling factors 
the energy per lattice site, the two-point and density-density correlations, and expectation values of 
powers of the on-site number operator determining the local atom number fluctuations (variance, 
skewness, kurtosis). We compare these expansions to numerical simulations of the infinite-size 
system to determine their range of applicability. We also discuss a new sum rule for the density- 
density correlations that can be used in both equilibrium and non-equilibrium systems. 


I. INTRODUCTION 


The Bose-Hubbard models capture key properties of numerous experimentally-relevant configurations of cold bosonic 
atoms placed in optical lattices [1-4]. The simplest of them is defined by the Hamiltonian 

H = ~J^2 (a\ +1 ai + h.c.) + \ X ”*(”* “ 1) > 

i i ( 1 ) 

= Sij, [ &i,aj ] =0, hi = fit a*, 

where the first term describes tunnelling between adjacent sites, while the second one accounts for on-site interactions. 
The competition between these two terms leads to the Mott insulator-superfluid quantum phase transition when the 
filling factor (the mean number of atoms per lattice site) is integer [5, 6]. The system is in the superfluid phase when 
the tunnelling term dominates ( J > J c ) whereas it is in the Mott insulator phase when the interaction term wins out 
(J < J c ). The location of the critical point depends on the filling factor n and the dimensionality of the system. We 
consider the one-dimensional model (1), where it was estimated that 


0.3 

for n = 1 


0.18 

for n = 2 . 

(2) 

0.12 

for n = 3 



It should be mentioned that there is a few percent disagreement between different numerical computations of the 
position of the critical point (see Sec. 8.1 of Ref. [4] for an exhaustive discussion of this topic). That affects neither 
our results nor the discussion of our findings. 

The Bose-Hubbard model (1), unlike some one-dimensional spin and cold atom systems [6, 7], is not exactly solvable. 
Therefore, it is not surprising that accurate analytical results describing its properties are scarce. To the best of our 
knowledge, the only systematic way of obtaining them is provided by the perturbative expansions [8-14]. In addition 
to delivering (free of finite-size effects) insights into physics of the Bose-Hubbard model, these expansions can be used 
to benchmark approximate approaches (see e.g. Refs. [15, 16]). 

We compute the following ground-state expectation values: the energy per lattice site E, the two-point correlations 
C(r) = (OjCLj+r), the density-density correlations D(r) = (hjfij+r), and the powers of the on-site number operator 
Q(r) = (h\) - n r . 

Our perturbative expansions are obtained with the technique described in Ref. [10] (see also Ref. [11] for a similar 
approach yielding the same results). The differences with respect to Ref. [10] are the following. First, we have 
computed perturbative expansions for the filling factors n = 2 and 3, which were not studied in Ref. [10]. Second, we 
have enlarged the order of all the expansions for the n = 1 filling factor that were reported earlier. Moreover, several 
perturbative results for the n = 1 case, that were not listed in Ref. [10], are provided in Appendix B. Third, we have 
computed perturbative expansions for the expectation values of different powers of the on-site atom number operator, 
which were not discussed in Ref. [10]. This allowed us for computation of the skewness and kurtosis characterizing 
on-site atom number distribution. Fourth, we have derived an important sum rule for the density-density correlations 
allowing for verification of all our perturbative expansions for these correlations. 

The range of validity of our perturbative expansions is carefully established through numerical simulations. There 
is another crucial difference here with respect to our former work [10]. Namely, instead of considering a 40-site system, 
we study an infinite system using the translationally invariant version of the Time Evolving Block Decimation (TEBD) 
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algorithm sometimes referred to as iTEBD [17] (where i stands for infinite). The ground state of the system is found 
by imaginary time propagation [18]. For the detailed description of the method and its relation to the density matrix 
renormalization group studies see the excellent review [19]. The application of iTEBD allows for obtaining results 
free of the finite-size effects from numerical computations (see Appendix A for the details of these simulations). Our 
symbolic perturbative expansions have been done on a 256 Gb computer. The numerical computations require two 
orders of magnitude smaller computer memory. 

The outline of this paper is the following. We discuss in Sec. II various identities that can be used to check 
the validity of our perturbative expansions. In particular, we derive there a sum rule for density-density correlation 
functions. Sec. Ill is focused on the ground state energy per lattice site. Sec. IV shows our results for the variance 
of the on-site atom number operator. Sec. V discusses expectation value of different powers of the on-site number 
operator and the related observables: the skewness and kurtosis of the local atom number distribution. Sec. VI 
discusses the two-point correlation functions. Sec. VII provides results on the density-density correlations. The 
perturbative expansions presented in Secs. Ill-VII are compared to numerics, which allows for establishing the range 
of their applicability. Additional perturbative expansions are listed in Appendices B, C, D for the filling factors 
n = 1,2, 3, respectively. The paper ends with a brief summary (Sec. VIII). 


II. GROUND STATE IDENTITIES AND SUM RULE 


There are several identities rigorously verifying our perturbative results. First, straight from the eigen-equation one 
gets that the ground state energy per lattice site, E, satisfies 

E = ~2JC(1) + D ^~ n _ (3) 

It is easy to check that our perturbative expansions - (8), (12), and (27) for n = 1; (9), (13), and (30) for n = 2; and 
(10), (14), and (33) for n = 3 - satisfy this identity. 

Combining this result with the Feynman-Hellmann theorem, 

dJ h \dJ /' 

we get 

p m = i jd. cm . 

A similar identity can be found in Sec. 7.1 of Ref. [4|. Once again, it is straightforward to check that our expansions 
for n = 1,2,3 satisfy this identity. 

Finally, we obtain a sum rule for the density-density correlations in a one-dimensional system 

OO 

D{ 0) - n 2 + 2 £ [ D ( r ) - n2 ] = ( 4 ) 


It is again an easy exercise to check that our expansions - (36)-(38) and (B6)-(B10) for n = 1; (39)-(41) and (C5)- 
(C7) for n = 2; and (42)-(44) and (D4)-(D6) for n = 3 - satisfy this sum rule 1 . Eq. (4) can be obtained from the 
sum rule for the zeroth moment of the dynamic structure factor (see Ref. [20] for a general introduction to a dynamic 
structure factor and its sum rules and Ref. [21] for their discussion in a Bose-Hubbard model). We have, however, 
derived it in the following elementary way. 

Consider a system of N atoms placed in the M -site periodic lattice ( N,M < oo). Assuming that the system is 
prepared in an eigenstate of the number operator, say | T), we have 


N 2 



f M 

\i= 1 


2 



M 


7 = 1 


(5) 


1 There is no need to perform the sum over infinite number of D(r)’s to see that our results satisfy the sum rule (4). This follows from 
the observation that D{r > 0 ) — n 2 = O (J 2r ). Thus, if our expansions for n = 1 (n = 2 and 3) are done up to the order J 16 (J 12 ), we 
need to know D(r) only for r = 1, 2, • ■ • ,8 (r = 1, 2, • • • ,6). 
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Figure 1: The energy per lattice site for different filling factors. Lines come from expansions (8)— (10), while dots show numerical 
results obtained using iTEBD code with the imaginary time evolution. Both here and in other figures we have (i) added blue 
dotted lines connecting the dots to facilitate quantification of the discrepancies between perturbative expansions and numerics; 
(ii) drawn red vertical dotted lines at the positions of the critical points; and (iii) used all the terms of the computed perturbative 
expansions listed in the paper to plot the perturbative results. 


The next step is to assume that the correlations ('&\hihj\'$ 1 } depend only on the distance between the two lattice sites. 
This assumption allows for rewriting Eq. (5) to the form 


D( 0) 


LM/2J 


r= 1 


D(r) 



= 0 , 


( 6 ) 


where [_;rj stands for the largest integer not greater than x, [M/2\ is the largest distance between two lattice sites in 
the M-site periodic lattice, and the prime in the sum indicates that in even-sized systems the summand for r = |_M/2J 
has to be multiplied by a factor 1/2. One obtains Eq. (4) by taking the limit of N, M —> oo such that the filling factor 
n = N/M is kept constant. Such a procedure is meaningful as long as the correlations D(r) tend to n 2 sufficiently 
fast as r increases, which we assume. The extension of the above sum rule to two- and three-dimensional systems is 
straightforward, so we do not discuss it. 

Instead, we mention that the sum rule (6) can be also applied to non-equilibrium systems satisfying the assumptions 
used in its derivation. It can be used either to study constraints on the dynamics of the density-density correlations or 
to verify the accuracy of numerical computations. Both applications are relevant for the studies of quench dynamics of 
the Bose-Hubbard model triggered by the time-variation of the tunnelling coupling J [22-24]. We mention in passing 
that a completely different work on the sum rules applicable to the Bose-Hubbard model can be found in Ref. [251. 

Finally, we mention that it has been shown in Ref. [10] that the ground state energy per lattice site and the 
density-density correlations in the Bose-Hubbard model are unchanged by the 


(7) 

transformation, while the two-point correlations transform under (7) as C(r) (—1 ) r C(r). Using the same reasoning 
one can show that Q(r ) is symmetric with respect to (7) as well. One can immediately check that all the expansions that 
we provide satisfy these rules. This observation provides one more consistency check of our perturbative expansions. 
Moreover, it allows us to skip the 0(J m+2 ) term by the end of every expansion ending with a J m term. 
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III. GROUND STATE ENERGY 


The ground state energy per lattice site E for the unit filling factor is 


E = 68 1267 U171 

4 9 81 1458 

32507578587517774813 t16 

466647713280000 ’ 


4902596 t12 

- J 

6561 


8020902135607 14 
2645395200 


while for n = 2 it is given by 


( 8 ) 


El t2 t4 49604 , 3385322797 

— - 3 J 2 + 8J 4 H-J 6 -, 

4 4 315 13891500 

7350064303936751836656911 12 
15282461406452625000 ’ 

and finally for n = 3 it reads 

* =5 - 67 2 + 31 J 4 + - 11^07105017 

44 63 36117900 

39433892936615327274896871074109109 12 

1229047086250770739427475000 ' 

The ground state energy for an arbitrary integer filling factor was perturbatively calculated up to the J 4 terms in 
Sec. 7.1 of Ref. [4]. Our expansions, of course, match this result. 

A quick inspection of Fig. 1 reveals that there is an excellent agreement between numerics and finite-order per¬ 
turbative expansions (8)-(10) not only in the whole Mott insulator phase, but also on the superfluid side near the 
critical point (see Ref. [16] for the same observation in the n = 1 system). This is a bit surprising for two reasons. 

First, it is expected that the perturbative expansions break down at the critical point in the thermodynamically- 
large systems undergoing a quantum phase transition. This, however, does not mean that our finite-order expansions 
(8)-(10) cannot accurately approximate ground state energy per lattice site across the critical point. 

Second, we find it actually more surprising that despite the fact that our finite-order perturbative expansions for 
both (7(1) and D{ 0) depart from the numerics on the Mott side, their combination (3) works so well across the critical 
point. The two-point correlation function (7(1) is depicted in Figs. 8-10, while D{ 0) is given by var(h) + n 2 , where 
var(n) is plotted in Fig. 2. It would be good to understand whether this cancellation comes as a coincidence due to 
the finite-order of our perturbative expansions (8)-(10). 


8 8232891127289 t10 

+ 168469166250 3 


(9) 


76233225199535567419 

3516204203386875 


10 


( 10 ) 


IV. VARIANCE OF ON-SITE NUMBER OPERATOR 

The most basic insight into the local fluctuations of the number of atoms in the ground state is delivered by the 
variance of the on-site number operator 


var(h) = (n 2 ) — (hi) 2 = D( 0) — n 2 . 


( 11 ) 


This quantity is experimentally accessible due to the spectacular recent progress in the quantum gas microscopy [26] . 
We find that for the unit filling factor 


... ot2 _ t4 2720 Tfi 70952 
var(n) =8J 2 - 24J 4 - —J 6 H-J 8 


+ 


9 81 

32507578587517774813^. 


16 


3888730944000 


176684 Tin 431428448 t12 104271727762891 t14 

_j iu j_j _j i4 

81 6561 330674400 


( 12 ) 


for the Filing factor n = 2 


/-x ,„„ T 4 396832 , 6770645594 t8 

varin ) =24J 2 - 192J 4 -J 6 +-J 8 

v ; 63 496125 

7350064303936751836656911 t12 

+-- J , 


32931564509156 10 
9359398125 


( 13 ) 


173664334164234375 
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Figure 3: Expectation values of the powers of the on-site number operator (15) for the unit filling factor. Lines show expansions 
(16)-(18), while the dots show numerics. 


and finally for n = 3 


var{n) =48 J 2 - 744J 4 - 


2946560 


J 6 + 


22414210034 


63 1289925 

39433892936615327274896871074109109 
13966444161940576584403125 


J 8 + 


609865801596284539352 
390689355931875 ' 


10 


J 


12 


(14) 


The comparison between these perturbative expansions and numerics is presented in Fig. 2. We see there that our 
expansions accurately match numerics in most of the Mott phase and break down near the critical point. It might 
be worth to note that these on-site atom number fluctuations are nearly the same at the critical point (2) for the 
different filling factors (they equal roughly 0.4 there). 
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Figure 4: Expectation values of the powers of the on-site number operator (15) for the n = 2 filling factor. Lines show 
expansions (19)—(21) , while the dots show numerics. 


V. POWERS OF NUMBER OPERATOR 


Further characterization of the fluctuations of the occupation of individual lattice sites comes from the study of 
expectation values of the integer powers of the on-site number operator 


Q(r) = {n\) - n r 


(15) 


for r > 2 (the r = 2 case was analyzed in Sec. IV). Once again, we mention that these observables can be experimentally 
studied [26]. 

For the unit filling factor, we get 


Q(3) =24J 2 - 56J 4 - 976J 6 + + + J 14 


81 

5938172375134531873121 . 
181474110720000 J 


243 


52488 


9258883200 


16 


(16) 


<3(4) =56 J 2 - 72 J 4 - 


22784 


-J b + 


355192 


81 




31533614 


729 


J w + 


16939285963 
26244 


J 12 + 


488931794121599 
661348800 


-J 


14 


12234501340429656667403 
116661928320000 


-J 


16 


(17) 


<3(5) =120 J 2 + 40 J 4 — 


18800 


-r + - 


601000 


81 


-r - 


123485195 


729 


J iU + 


31523026139 


17496 


- J - 


1978940191363981 
1322697600 ' 


•14 


2143214705361163325357 


6666395904000 
For two atoms per site, we obtain 


-J 


16 


<3(3) =144J 2 - 1072 J 4 - J* + 1770730207 _ 436 J 8 - 077lll,:W - ,ill,i, ^' 7 ' . / :-, 


49 


24310125 


2162020966875 


15451331550936239672643340032833 
60371453478515288484375 


-J 1 


(18) 


(19) 


cnn /2 or70Cj4 1885928848 76 , 5417457952036 t8 107844070676948560562 t 
CJ (4) —600J 3736 J -- txtt— <7 —t—-.... _ J - , —, _....... _- J 


10 
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40516875 
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J 
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(20) 
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Figure 5: Expectation values of the powers of the on-site number operator (15) for the n — 3 filling factor. Lines show 
expansions (22)-(24), while the dots show numerics. 


Q{ 5) =2160 J 2 - 9440J 4 - 


735 24310125 

480125387136897585787036245853433 12 

+ 120742906957030576968750 ' ' 


Finally, for three atoms per site we derive 


Q{ 3) =432 J 2 - 6472 J 4 - 


1284712 Tfi 1195336576618 tS 
+ 16769025 + 


3 16769025 

1273450413079818438111858514006273177409357 
50089814000150161485366743625000 ' 


zrozuzuyooo/o 


27678339796268712326815412 10 


3184508940200713125 


12 


( 21 ) 


( 22 ) 


Q{ 4) =2640J 2 - 36400J 4 - 


846750928 


-J b - 


59064210154568 


r - 


315 23476635 

66279835521862060615675760372212019355789667 12 

425763419001276372625617320812500 ’ 


1031160890254623471701872 10 
974849675571646875 


(23) 


Q(5) =13680J 2 - 163280J 4 - 




J b - 


i 23476635 

9189183527664354899691980380144063394455799 12 
11353691173367369936683128555000 


9553526820602139375 


(24) 


These expansions are compared to numerics in Figs. 3-5. They reproduce the numerics in the Mott insulator phase 
in the same range of the tunneling coupling J as our expansions for the variance of the on-site number operator. 

Using expansions (16)-(24) one can easily go further, i.e., beyond the variance, in characterization of the on-site 
atom number distribution. For example, one can easily compute the skewness [27, 28] 


= {{hi - n) 3 ) 
<0 hi - n) 2 ) 3 / 2 

and the kurtosis [27, 28] (also referred to as excess kurtosis) 

{{hi ~ nf) 
{{hi - n) 2 ) 2 


(25) 


K = 


(26) 



























Figure 6: The skewness of the on-site atom number distribution. Lines show Eq. (25) computed with expansions from Secs. 
IV and V. Dots show numerics. 


The skewness is a measure of a symmetry of the distribution. It is zero for a distribution that is symmetric around 
the mean. We plot the skewness in Fig. 6 and find it to be positive in the Mott insulator phase, which indicates that 
the distribution of different numbers of atoms is tilted towards larger-than-mean on-site occupation numbers. This 
is a somewhat expected result given the fact that the possible atom occupation numbers are bounded from below by 
zero and unbounded from above. Given the fact that 151 < 1/2 in Fig. 6, one may conclude that the on-site atom 
number distribution is “fairly symmetric” in the Mott phase according to the criteria from Ref. [28]. 

The kurtosis quantifies whether the distribution is peaked or flat relative to the normal (Gaussian) distribution. 
It is calibrated such that it equals zero for the normal distribution of arbitrary mean and variance. K > 0 ( K < 0) 
indicates that the studied distribution is peaked (flattened) relative to the normal distribution. We plot the kurtosis in 
Fig. 7. As J —> 0 one easily finds from our expansions that K ~ J~ 2 . This singularity reflects the strong suppression 
of the local atom number fluctuations in the deep Mott insulator limit. The kurtosis monotonically decays in the 
Mott phase (Fig. 7). 

To put these results in context, we compare them to the on-site atom number distribution in the deep superfluid 
limit of J —> oo (the Poisson distribution [29]). The probability of finding s atoms in a lattice site is then given in the 
thermodynamic limit by exp(— n)n s /s\ : where n is the mean occupation. One then finds that S = 1 /y/n and K = 1/n 
for the Poisson distribution. Keeping in mind that the Gaussian distribution is characterized by S = K = 0, we can 
try to see whether the on-site atom number distribution near the critical point is Gausssian-like or Poissonian-like. 

We see from Figs. 6 and 7 that at the critical point (2) we have S ss 0.22,0.11,0.07 and K « 0.19,0.3,0.4 for 
n = 1,2,3, respectively. Therefore, the real distribution lies somehow between Poissonian and Gaussian. The skewness 
suggests that for these filling factors the distribution at the critical point is more Gaussian than Poissonian. On the 
other hand, the kurtosis for n = 1 (n = 2,3) is more Gaussian (Poissonian). From this we conclude that for the 
unit filling factor the on-site atom number distribution at the critical point is better approximated by the Gaussian 
distribution. 


VI. TWO-POINT CORRELATIONS 

The two-point correlation functions play a special role in the cold atom realizations of the Bose-Hubbard model 
[30-32]. Their Fourier transform provides the quasi-momentum distribution of a cold atom cloud, which is visible 
through the time-of-flight images that are taken after releasing the cloud from the trap. 
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Figure 7: The kurtosis of the on-site atom number distribution. Lines show Eq. (26) computed with expansions from Secs. IV 
and V. Dots show numerics. 
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Figure 8: The two-point correlation functions for the unit filling factor. Lines from top to bottom correspond to r = 1, 2,3, 
respectively. They depict perturbative expansions (27)-(29). The numerics is presented with dots. 


For the filling factor n = 1, they are given by 
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Figure 9: The two-point correlation functions for the n — 2 filling factor. Lines from top to bottom correspond to r = 1, 2,3, 
respectively. They depict perturbative expansions (30)-(32). The numerics is presented with dots. 
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Figure 10: The two-point correlation functions for the n — 3 filling factor. Lines from top to bottom correspond to r = 1, 2,3, 
respectively. They depict perturbative expansions (33)-(35). The numerics is presented with dots. 


no r3 8324 t 5 126040 t7 7883333 t9 220980576341 

C( 3) =88 J 6 -J 5 -I-J 7 -I-- 

w 9 81 486 1049760 

82283484127688477 13 
+ 61725888000 ’ 


and for the filling factor n = 2 they are 
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and finally for n = 3 they read 
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(34) 


(7(3) =2928J 3 
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(35) 


Expansions up to the order J 3 for (7(1), (7(2), and (7(3) at arbitrary integer filling factors are listed in Sec. 7.1 of 
Ref. [4] and agree with our results. 

We see in Figs. 8-10 that the above perturbative expansions break down within the Mott insulator phase (the 
larger r, the deeper in the Mott phase the expansion breaks down). We notice that it is instructive to compare the 
value of the correlations <7(r) around the critical point to their deep superfluid limit. C(r) in the J —> oo limit tends 
to n (see e.g. Appendix B of Ref. [10]). Therefore, the three correlation functions C(r = 1, 2, 3) reach at least 50% of 
their deep superfluid value near the critical point, which well illustrates the significance of quantum fluctuations at 
the critical point. 

The ground state quasi-momentum distribution is defined as 


n(k) 


1 

M 


M 

{a^cLs) exp [ik(m 

m,s= 1 


*)]■ 


where M stands for the number of lattice sites (we skip the prefactor proportional to the squared modulus of the 
Fourier transform of the Wannier functions; see Ref. [31] for details). Taking the limit of M —> oo at the fixed integer 
filling factor n, one gets 


OO 

n(k) = n + 2 ^ (7(r) cos (rk). 

r —1 

Using Eqs. (27)-(29) and (B1)-(B5) for n = 1, Eqs. (30)-(32) and (C1)-(C4) for n = 2, and Eqs. (33)-(35) and 
(D1)-(D3) for n = 3 the state-of-the-art high-order perturbative quasi-momentum distributions for different filling 
factors can be obtained. These results can be compared to Ref. [14], where an expression with terms up to J 3 for an 
arbitrary filling factor is computed. As expected, we find these results in agreement with our findings. 


VII. DENSITY-DENSITY CORRELATIONS 


Similarly as the observables from Secs. IV and V, the density-density correlations can be experimentally approached 
through the technique discussed in Ref. [26]. 

The density-density correlations are given for n = 1 by 
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Figure 11: The density-density correlation functions for the unit filling factor. Lines from bottom to top illustrate perturbative 
results for r = 1,2,3, respectively. Dots show numerics. Perturbative expansions are given by Eqs. (36)-(38). 



J 

Figure 12: The density-density correlation functions for the n = 2 filling factor. Lines from bottom to top illustrate perturbative 
results for r = 1,2,3, respectively. Dots show numerics. Perturbative expansions are given by Eqs. (39)-(41). 
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1588041877 T , n 1710030328933 , 2 
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(38) 


for n = 2 they read 


11(1) =4 — 12 J 2 + 
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J 4 — 


80553632 

33075 


J b - 
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121550625 


J® - 
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6486062900625 


J 


10 


+ 


331686439652436848222471319678887 
72445744174218346181250000 


J 


12 


(39) 


812 4 189419192 6 140772979852859 8 569733162420769673609 jl0 

~~ J + 11025 364651875 + 324303145031250 


8056033392249986400009146407648 

503095445654294070703125 


12 


(40) 
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Figure 13: The density-density correlation functions for the n = 3 filling factor. Lines from bottom to top illustrate perturbative 
results for r = 1, 2, 3, respectively. Dots show numerics. Perturbative expansions are given by Eqs. (42)-(44). 
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2028480836878113693075000000 
and finally for n = 3 they can be written as 


-J 1 
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71(1) —9 — 24J H-- —J — _ —J —- ^ ^ ^ _ ■' —_ ____ - J 


6615 


4108411125 


200624063232644926875 


1542480719505910230376228320731168033769541 

3193225642509572794692129906093750 


(42) 


J 1 


3160 4 294843464 6 76035818562996449 8 83670951564711862884744588592 10 

3 + 2205 12325233375 + 1003120316163224634375 

143706393091669463828236051561683582721397 12 

132705481247151077182010593500000 ’ 


(43) 


12148432 6 246576902129764 8 1185488040768577918685665169 10 

135 + 14586075 926242212523753125 

13080624640958701853202057691706510935349239 12 

+ 285200376364491350084145573750000 ' ' 


(44) 


The correlation functions D{ 1) and D(2) were computed for an arbitrary integer filling factor up to the order J 4 in 
Sec. 7.1 of Ref. [4], These results agree with our expansions. 

The comparison between our perturbative expansions and numerics is presented in Figs. 11-13 for different filling 
factors. The expansions break down near the critical point on the Mott side of the transition. Comparing Figs. 8-10 to 
Figs. 11-13, we see that expansions for the two-point and density-density correlations break down in similar distances 
from the critical point. Moreover, this comparison shows that the two-point correlations change more appreciably 
within the Mott phase than the density-density correlations. We attribute it to the constraints that are imposed on 
the density-density correlations due to the atom number conservation. 


VIII. SUMMARY 

We have computed state-of-the-art high-order perturbative expansions for several observables characterizing ground 
state properties of the one-dimensional Bose-Hubbard model in the Mott phase. As compared to our former results 
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for the filling factor n = 1 [10], we have extended our analysis by considering the filling factors n = 2 and 3 (we have 
also enlarged the number of terms for the n = 1 case). We have characterized the on-site atom number distribution 
by giving the predictions for the skewness and kurtosis. Those may serve as useful benchmarks for experimental 
in-situ observations [26] . We have also derived in a simple way an important sum rule applicable to both equilibrium 
and non-equilibrium density-density correlations. That sum rule allows for verification of our perturbative expansions 
and it may be useful for checking the consistency of experimental data. We have also carefully established the 
range of applicability of our perturbative expansions through numerical simulations. The expansions discussed in 
this work can be easily typed or imported into computer software such as Mathematica or Maple and used for 
benchmarking approximate approaches, comparing theoretical predictions to experimental measurements, testing 
Pade approximations, etc. 
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Appendix A: iTEBD simulations 

There are two factors that have to be taken care of to assure the convergence of results in the numerical implemen¬ 
tation of iTEBD. The first one is the maximal allowed number of bosons per site assumed in the variational ansatz, 
Nmax • We take N max = 6 for the filling factor n = 1 up to N max = 12 for n = 3. We have checked that these values 
lead to converged results. The second important factor is the number of Schmidt decomposition eigenvalues, x, kept 
during each step of the procedure [17, 19]. x ma y be quite small deep in the Mott regime (of about 20) while it must 
be significantly increased close to the critical point and in the superfluid regime. We have found that for reliable 
energy, particle number variance, as well as two-point correlations with small r the choice of x = 150 was largely 
enough (with the relative error of the order of 10 - ' in energy and 10 -5 in particle number variance). Let us note that 
the numerical studies of long-range correlations (r of the order of a hundred) require taking x > r at least [18]. 


Appendix B: One atom per site 


Our remaining perturbative expansions for the n = 1 filling factor are listed below. 
The two-point correlations: 


(7(4) = 450J 4 


186608 


27 


-J b + 


7565704 

243 


J 8 + 


1493509507 10 
17496 


858313783040137 12 
440899200 


(Bl) 


(7(5) = 2364J 5 


3894512 j7 


250517014 q 

- J 

729 


25842700043 n 
209952 


(B2) 


(7(6) = 12642J 6 


78008768 
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J 8 + 


6836492080 10 
2187 


(B3) 


(7(7) = 68464J 7 


1522020908 

729 


(B4) 


(7(8) = 374274J 8 . 


The density-density correlations: 


n 1 741706 r8 ( 93328235 Tl0 288653212433561 rl2 
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-16 
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1738824899 10 91423339623697 12 897923823504590743589 14 

1 8748 + 8817984 3888730944000 

40404128939318210395355039327 16 

+ 15091387047475200000 ’ 


(B7) 


D( 6) =1 


369347437555 12 22768945554355275259 14 

78732 + 77774618880 


20246612891148030348297322711 16 
2515231174579200000 


(B8) 


£>( 7 ) =1 


1771595060952703 14 
15116544 


4869453765809764188858679 16 
571643448768000 


(B9) 


D( 8) =1 


415126490285461535 16 
136048896 


(BIO) 


Appendix C: Two atoms per site 


Our remaining perturbative expansions for the n = 2 filling factor are listed below. 
The two-point correlations: 


C(4) = 6450J 4 


7683200 Tfi 291093977979158 t8 10163319115920107956583 Tin 

_/ b j_/ 8 j_7 1U 
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The density-density correlations: 
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D(6) =4- 
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Appendix D: Three atoms per site 


Our remaining perturbative expansions for the n = 3 filling factor are listed below. 
The two-point correlations: 


C(4) = 35700J 4 
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The density-density correlations: 
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